function qs = ElnpY_CZmu_losses(mix, data, vars)

% E[ln p(Y|C, Z, mu)]

K = mix.ncentres;
[N D] = size(data.Y);
Dlog2pi = D*log(2*pi);
m = mix.varposterior.m;
T = mix.varposterior.T;

% precomputed values and vectorised data
invC_vec = reshape(data.invC, D*D, N);
logdetC = data.logdetC;
invCy = data.invCy;
yinvCy = data.yinvCy;

% vectorised model params
mmt = zeros(D, D, K);
for k = 1:K
    mmt(:, :, k) = m(k, :)'*m(k, :);
end
mmt_vec = reshape(mmt, D*D, K);
T_vec = reshape(T, D*D, K);

% compute the expectation using vectorised form
q = -Dlog2pi - logdetC - yinvCy;    % yinvCy is part of the quadratic form
q = q*ones(1, K);
q = q - invC_vec'*mmt_vec + 2*invCy*m' - invC_vec'*T_vec;

qs = zeros(1, K);
for k = 1:K
    comps = [1:(k-1) (k+1):K];
    Z = vars.Z(:, comps) + realmin;
    Z = Z./(sum(Z, 2)*ones(1, K-1));
    
    qs(k) = 0.5*sum(sum(Z.*q(:, comps)));
end
